#------------------------------------------------------------------------------

#  读取消费过的会员信息：rawdata
data<-rawdata

# 定义分类变量和数值变量
factor.names<-names(data)[!(sapply(data,class) %in% c('numeric','interger'))]
vector.names<-names(data)[sapply(data,class) %in% c('numeric','interger')]

# 数据初步整理

data$memberid<-data$MemberID<-data$Birthday<-data$email<-NULL

data<-data[outlier.delete(data[,vector.names]),]

# 缺失值处理
data<-lossdelete(data)
lossrate(data)

# add.miss 自定义插补数值型变量
num
data<-add.miss(data,names(data)[lossrate(data[,vector.names])>0],num)

# 无法插补的缺失值去掉
data<-na.omit(data)

# 相关分析
library(psych)
mem.vcor<-corr.test(data[,vector.names],use="pairwise")

# 主成份数据集pca.data，用以评价推荐指数
pca.data<-data[,vector.names[grep("mgm",vector.names)]]

# 看指标之间的相关性
data.cor<-cor(pca.data,pca.data,"pairwise.complete.obs")

# 剔除与其他指标间相关性小的指标
vname<-names(pca.data)[apply(data.cor,2,function(i) max(abs(i[i<1]))>0.1)]

#再看剔除指标后的相关性
data.cor<-cor(data[,vname],data[,vname],"pairwise.complete.obs")        

pca.data<-scale(pca.data[,vname])

mem.pca<-princomp(pca.data,cor=T,scores=T)
summary(mem.pca,loadings=F)

#画主成分的碎石图
screeplot(mem.pca,type =  "lines")

#得到主成分的值
pca.score<-predict(mem.pca)[,mem.pca$sdev>1]
pca.score<-((mem.pca$sdev^2)/sum(mem.pca$sdev^2))[mem.pca$sdev>1] %*% t(pca.score)
pca.score<-as.vector(t(pca.score))
summary(pca.score)

data$score<-pca.score

data$response<-cut(pca.score,breaks=c(min(pca.score),-1,0,1,max(pca.score)+0.1),labels=c(-2,-1,1,2))
table(data$response)
prop.table(table(data$response))

pca.mean<-NULL
for(i in vector.names[grep("mgm",vector.names)]){
  pca.mean<-rbind(pca.mean,tapply(data[,i],data$response,mean))
}
rownames(pca.mean)<-vector.names[grep("mgm",vector.names)]
pca.mean<-round(pca.mean,6)
pca.mean


# 建训练样本
sampleid<-sample(1:nrow(data),round(nrow(data)*0.8,0))
traindata<-data[sampleid,!(names(data) %in% colnames(pca.data))]
testdata<-data[-sampleid,]
traindata$score<-NULL

# 建平衡样本
num<-round(min(table(traindata$response)$freq),0)
train<-data.balance(traindata,"response",c(-2,-1,1,2),c(1.2,2.5,2.5,1.2)*num)

# 模型训练
library(randomForest)
formula<-as.formula("response~.")

model<-train.model(traindata,formula,20,testdata)

n<-which(sapply(model,function(x) x$value[2])==max(sapply(model,function(x) x$value[2])))[1]

mem.rf<-model[[n]]$model

varImpPlot(mem.rf)
importance(mem.rf)


model<-train.model(train,formula,20,testdata)

n<-which(sapply(model,function(x) x$value[2])==max(sapply(model,function(x) x$value[2])))[1]

mem.rf2<-model[[n]]$model

varImpPlot(mem.rf2)
importance(mem.rf2)


# 模型评价

modelA<-model.evaluation(data,response,mem.rf)
modelB<-model.evaluation(data,response,mem.rf2)

mem.rf<-iif(modelA[4]>model[4],mem.rf,mem.rf2)

rm(list=ls()[!(ls() %in% c("mem.rf","rawdata"))])























#==================================================================================
#==================================================================================
#   Appendix (Functions)
#==================================================================================
#==================================================================================

#---查看缺失率,删除缺失率高的字段
lossrate<-function(x){
  a<-apply(is.na(x),2,mean)
  a<-round(a,3)
  return(a)
}

lossdelete<-function(x){
  a<-apply(is.na(x),2,mean)
  a<-round(a,3)
  deletecols<-names(a[a>0.8])
  deleterows<-names(a[a>0&a<0.2])
  for(i in deletecols){
    x[,i]<-NULL
  }
  for(j in deleterows){
    x<-x[!is.na(x[,j]),]  
  }
  return(x)
}

#---分类型变量数据集，相关性分析
cor.factor<-function(data){
  n<-names(data)
  cor<-NULL
  p<-NULL
  for(i in n){
    for(j in n){
      newdata<-na.omit(data[,c(i,j)])
      chisq<-chisq.test(table(newdata[,1],newdata[,2]))
      cor<-c(cor,chisq[[1]])
      p<-c(p,chisq[[3]]>0.05)
    }
  }
  cor[p]<-0
  cor<-matrix(cor,nrow=length(n))
  colnames(cor)<-n
  rownames(cor)<-n
  return(cor)
}


#---删除没有意义的字段
valid.data<-function(data){
  invalid<-na.omit(names(data)[sapply(data,sd,na.rm=T)==0])
  valid<-setdiff(names(data),invalid)
  data<-data[,valid]
  return(data)
}


#---插补缺失值
add.miss<-function(data,col,num){
  if(length(col)==1){
    i=col
    data[is.na(data[,i]),i]<-num
  }
  else{
    for(i in col){
      data[is.na(data[,i]),i]<-num[i]
    }
  }
  return(data)
}

#---oneway 方差分析
fc.test<-function(factor,varible,data){
  fc=NULL
  for(i in factor){
    for(j in varible){
      formula=as.formula(paste(j,"~",i))
      fc=c(fc,oneway.test(formula,data=data,var.equal = TRUE)$p.value<0.05)
    }
  }
  fc=matrix(fc,nrow=length(varible),ncol=length(factor),byrow=F,dimnames=list(varible,factor))
  return(fc)
}

#---创建平衡样本
data.balance<-function(data,v,vlevels,num){
  id=NULL
  for(i in vlevels){
    N<-which(data[,v]==i)
    nsample<-sample(c(1:length(N)),num[i],replace=T)
    id<-c(id,N[nsample])
  }
  newdata<-data[id,]
  return(newdata)
}

#---分类型变量数据集，相关性分析
cor.factor<-function(data){
  n<-names(data)
  cor<-NULL
  p<-NULL
  for(i in n){
    for(j in n){
      newdata<-na.omit(data[,c(i,j)])
      chisq<-chisq.test(table(newdata[,1],newdata[,2]))
      cor<-c(cor,chisq[[1]])
      p<-c(p,chisq[[3]]>0.05)
    }
  }
  cor[p]<-0
  cor<-matrix(cor,nrow=length(n))
  colnames(cor)<-n
  rownames(cor)<-n
  return(cor)
}

#---异常值处理

outlier.delete<-function(data){
  ks<-sapply(data,function(x) ks.test(x+runif(length(x),-0.01,0.01),"pnorm",mean(x),sd(x))$p.value)
  mean<-sapply(data,mean,na.rm=T)
  sd<-sapply(data,sd,na.rm=T)
  Q1<-sapply(data,quantile,0.25,na.rm=T)
  Q3<-sapply(data,quantile,0.75,na.rm=T)
  QR<-Q3-Q1
  
  outlier<-apply(data[,ks<0.05],1,function(x) mean(x>Q1[ks<0.05]-3*QR[ks<0.05] & x<Q3[ks<0.05]+3*QR[ks<0.05],na.rm=T))
  newID<-rownames(data)[outlier>quantile(outlier,0.2)]
  
  if (sum(ks>0.05)>0){
    outlier<-apply(data[,ks>0.05],1,function(x) mean(x>mean[ks>0.05]-3*sd[ks>0.05] & x<mean[ks>0.05]+3*sd[ks>0.05],na.rm=T))  
    newID<-rownames(data)[outlier>quantile(outlier,0.2)]
  }
  return(newID)  
}


#---模型训练

train.model<-function(traindata,formula,times,testdata){
  model<-list()
  for(n in 1:times){
    mem.rf<-randomForest(formula,data=traindata,na.action=na.omit,ntree=120)
    testpre<-predict(mem.rf,testdata)
    testpretable<-table(testdata$response,testpre)
    accuracy<-length(which(testdata$response==testpre))/nrow(testdata)
    target<-sum(testpretable['2','2']+testpretable3['-2','-2'])/sum(testpretable3['2',]+testpretable3['-2',])*sum(testpretable3['2','2']+testpretable3['-2','-2'])/sum(testpretable3[,'2']+testpretable3[,'-2'])
    model[[n]]<-list(value=c(accuracy,target),model=mem.rf)  
  }
}


#---模型评价

model.evaluation<-function(data,response,model){
  testpre<-predict(model,data)
  testpretable<-table(data$response,testpre)
  accuracy<-length(which(data$response==testpre))/nrow(data)
  target.accuracy<-sum(testpretable['2','2']+testpretable['-2','-2'])/sum(testpretable['2',]+testpretable['-2',])
  target.recall<-sum(testpretable['2','2']+testpretable['-2','-2'])/sum(testpretable[,'2']+testpretable[,'-2'])
  return(c(accuracy,target.accuracy,target.recall,accuracy*target.accuracy*target.recall))
}
